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Abstract 

In this article we describe a stable partitioned algorithm that overcomes the added mass instability arising in fluid- 
structure interactions of light rigid bodies and inviscid compressible flow. The new algorithm is stable even for 
bodies with zero mass and zero moments of inertia. The approach is based on a local characteristic projection of 
the force on the rigid body and is a natural extension of the recently developed algorithm for coupling compressible 
flow and deformable bodies [TJ El [3] . Normal mode analysis is used to prove the stability of the approximation for 
a one-dimensional model problem and numerical computations confirm these results. In multiple space dimensions 
the approach naturally reveals the form of the added mass tensors in the equations governing the motion of the rigid 
body. These tensors, which depend on certain surface integrals of the fluid impedance, couple the translational and 
angular velocities of the body. Numerical results in two space dimensions, based on the use of moving overlapping 
grids and adaptive mesh refinement, demonstrate the behavior and efficacy of the new scheme. These results include 
the simulation of the difficult problem of a shock impacting an ellipse of zero mass. 

Keywords: fluid-structure interaction, added mass instability, moving overlapping grids, compressible fluid 
flow, rigid bodies 
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1. Introduction 

An important class of fluid-structure interaction (FSI) problems are those that involve the interaction of moving 
bodies with high-speed compressible fluids. For example, understanding the impact of shock or detonation waves on 
rigid structures and embedded rigid bodies is of great interest. The numerical simulation of such problems can be 
difficult, and many techniques have been developed to address various facets of the problem. For a review of FSI 
see [4] for example. One particularly challenging aspect has been the presence of numerical instabilities that can arise 
when simulating problems with light bodies. This so-called added-mass instability is associated with the fact that the 
reaction of a body to an applied force depends not only on the mass of the body but also on the fluid displaced by the 
body through its motion. Traditional partitioned FSI schemes do not take into account the strong coupling between 
the fluid and solid and thus can exhibit an instability whereby the over-reaction of a light solid to an applied force 
from the fluid leads in turn to an even larger reaction from the fluid and so on. Fully coupled monolithic approaches 
to FSI can overcome the unstable behavior but are generally more expensive, can be difficult to implement, and 
may require advanced solvers or preconditioners. For compressible fluids the instability in partitioned algorithms 
can often be suppressed by choosing a smaller time-step (as the analysis in this article demonstrates). However, the 
stable time-step goes to zero as the mass of the body goes to zero and thus alternative approaches to removing the 
instability are desirable. 

In a recent series of articles, we have developed a set of stable interface approximations for partitioned solutions 
procedures that couple compressible fluids and deformable bodies [T] [2] |3] . In [T] [2] the interface approximation is 
based on a local characteristic analysis that results in an impedance weighted projection of the velocity and forces on 
the interface. These methods ensure the stability of the partitioned FSI scheme even for light solids. In this article 
we extend these ideas to the coupling of compressible fluids and rigid bodies. The key idea presented in this article 
can be introduced by considering the equations of motion for a rigid body (the full set of equations are presented in 



detail in Section 2.2| 



m&v 6 = T, (1) 
Alj = -uj x (Auj) + T, (2) 

where m& is the mass of the body, V(,(t) is the velocity of the center of mass, u(t) the angular velocity and A the 
moment of inertia tensor. T and T are, respectively, the force and torque on the body arising from the fluid forces 
on the surface of the body. From Equations |l])-(|2| it would at first seem impossible to solve for V;, and/or lj when 
nib = and/or A — 0, as the equations apparently become singular. However, from a local characteristic analysis of 
the appropriate fluid-structure Riemann problem, we can determine how T and T implicitly depend on the motion 
of the body, 

T = -A vv v b - A vlj uj +3F, T = -A" v v b - A uw u + f. (3) 

The matrices A 13 are the added-mass tensors; these are defined in terms of certain integrals of the fluid impedance 
over the boundary of the rigid body (see Section [6b. It is worth pointing out that the concept of added-mass has 
a long history in describing the motion of embedded bodies in both compressible and incompressible flows. For the 
compressible regime the recent article [5] nicely discusses the history as well as modern developments. 

Using the form of Equation |3| as a starting point, we define a partitioned FSI scheme that remains stable with 
a large time-step (i.e. the usual time-step restriction associated with the fluid domain in isolation) even as m b or A 
go to zero, provided the added-mass tensors satisfy certain properties. This approach relies on the use of an implicit 
time stepping method for the evolution of the rigid body. The new added-mass scheme is analyzed in detail for a one- 
dimensional model problem consisting of a rigid body embedded in a fluid governed by the linearized Euler equations. 
Both a first-order accurate upwind scheme and the second-order accurate Lax-WendrofF scheme are analyzed using 
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normal mode stability theory [6]. When the rigid body is integrated with an A-stable time-stepping method, the 
resulting partitioned FSI scheme is shown to be stable with a large time step even when m b — 0. 

The added-mass scheme is implemented in two space dimensions using the moving overlapping grid technique 
described in [7|. In this approach, local boundary fitted curvilinear grids are used to represent the bodies and these 
move through static background grids that are often chosen to be Cartesian grids for efficiency. Adaptive mesh 
refinement (AMR) is used on both curvilinear and Cartesian grids to dynamically increase resolution locally in space 
and time. We solve the compressible Euler equations on possibly moving grids in the fluid domain using a high- 
order extension of Goudnov's method. The Newton-Euler equations (with added-mass corrections) are solved for the 
motion of the rigid-body using an implicit Runge Kutta scheme (in contrast to the explicit time-stepping method 
used previously in [7]). 

In general, the added-mass scheme proposed here could be used in conjunction with any number of FSI approaches. 
The treatment of moving geometry is a major component for coupling fluid flow to the motion of rigid bodies and 
many techniques have been considered. One class of methods relies on a fixed underlying grid for the fluid domain and 
includes, embedded boundaries |5j, immersed boundaries [5][TU], level sets [111112] . and fictitious domain methods 13 . 
A second class of methods uses body conforming meshes and allows the mesh to deform in response to the motion of 
the body. Popular in this class of methods are ALE [1411151 [151117] . multiblock [TH], and general moving unstructured 
grids H5J. 

The remainder of this article is structured as follows. In Section[2j the governing equations of inviscid compressible 
flow for the fluid, and the Newton-Euler equations for rigid body motion are presented. Section [3] provides some 
motivation for, and the derivation of, our interface projection scheme in one dimension, showing the origin of the 
added-mass terms in the equation of motion for the rigid body. In Section[4]this approximation is incorporated into a 
partitioned FSI scheme for a one-dimensional FSI model problem. The stability of this new added-mass scheme, as well 
as the traditional coupling scheme, is analyzed using normal mode theory. Section [5] provides numerical confirmation 
of the theoretical results for the one-dimensional problem, demonstrating the expected convergence rates and stability 
properties. Extension of the algorithm to multiple space dimensions is presented in Section [6] showing the derivation of 
added-mass tensors. The time-stepping procedure for the overlapping grid FSI algorithm is summarized in Section[7| 
Results for two-dimensional problems are presented in Section [8] These include a smoothly receding rigid piston with 
known solution, a smoothly accelerated ellipse which is compared to the traditional algorithm and a shock-driven 
zero mass ellipse. The last example is also used to demonstrate the use of adaptive mesh refinement (AMR) and is 
particularly challenging and interesting. Concluding remarks are given in Section [9] In |Appcndix A we derive the 
exact solutions used in the numerical verification of the one-dimensiomal model problem. Finally in | Appendix B"] we 
present the form of the added mass matrices for a number of simple shapes in two and three dimensions. 



2. Rigid bodies and compressible flow in multiple space dimensions 

In this section we define the governing equations for the fluid domains and the rigid bodies. The equations are 
presented in three space dimensions which serves as a general model. Simplifications to one and two space dimensions, 
as well as linearization, will be performed later as appropriate. 

2.1. The Euler equations for an inviscid compressible fluid 

We consider the evolution of a compressible inviscid fluid with an embedded rigid body. The governing equations 
for the fluid domain !!cl 3 are the compressible Euler equations 

d t vf + V • f (w) = 0, xefi, t > 0, (4) 

where w = [p, pv,p£] T is the vector of conserved variables (density, momentum, energy), v is the velocity, and 
f = [pv, pv ® v + pi, (p£ + p)v] T is the flux. The total energy is given by p£ — p/(j — 1) + |pjv| 2 assuming an ideal 
gas with a constant ratio of specific heats. 

2.2. The Newton-Euler equations for the motion of a rigid body 

The equations of motion for the rigid body are the Newton-Euler equations which can be written as 



x 6 = Vb, (5) 

m b ir b = T, (6) 

Aw = -WAw + T, (7) 

E = WE. (8) 



Here m b is the mass of the body, Xb(t) G R 3 is the position of the center of mass, and Vb(t) G R 3 is the velocity of 
the center of mass. The moment of inertia matrix A £ R 3x3 is defined by 

A(t)= p 6 (x) y T y/-yy T dx, y = x-x 6 , 

JB{t) 1 J 
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where pb(x) defines the mass density of the body and B(t) C R 3 defines the region occupied by the body. The inertia 
matrix is symmetric and positive semi-definite (positive definite if p&(x) > 0) and can be written in terms of the 
orthogonal matrix E £ R 3x3 , whose columns are the principle axes of inertia, ej(t), and the diagonal matrix A whose 
diagonal entries are the moments of inertia, U, 

A = EAE T , £ , = [e 1 e 2 e 3 ], Ae* = ije,, A = diag(/i, h,I 3 ), ef e 3 - = <5y. 

The angular momentum of the body is h = Aijj where oj(t) £ R 3 is the angular velocity. The matrix W in Q is the 
angular velocity matrix given by 



W = Cross(w) = 





— OJ2 Wl 



-0J3, UJ2 
-Wl 





( i.e. Wa. = w X a) 



The total force and torque on the body are given by 



T 



I i s ds + fb, (f s = surface forces, ft,= body force), 

J dB 



(x — Xj) x f s ds + g;,, (torque), 



(9) 

(10) 
(11) 



dB 



Given J-(t) and T(t), along with initial conditions, X(,(0), Vb(0), t^(0), and -E(O), equations ([5|-((8| can be solved 
to determine X{,(t), ~Vb(t), <^(i), and E(t) as a function of time. 

The motion of a point r(i) attached to the body is given by a translation together with a rotation about the 
initial center of mass, 

r(i) = x b (t) + R(t)(r(0) -K b (Q)), 
where R(t) is the rotation matrix given by 



R(t) = E(t)E T {Q). 



(12) 



The velocity of this point is 



r(i) = Vi(t) + Wfl(t)(r(0) - x b (0)), 
= v i (t)+W(r(t)-x i (t)) ) 
= Vi,(t) +w x (r(t) -Xi,(t)), 

Letting y = y(r) = r(t) — Xb(t) it follows that the velocity of the point r can be written in the form 

r(t) = v b (t) - Yu, 

where Y(t) is the matrix 



Y = Cross(y) 






-2/3 


2/2 


2/3 





-2/1 


-V2 


2/1 






(13) 



(14) 



2. 3. The coupling conditions for rigid bodies and inviscid compressible flow 

On an interface between a fluid and a solid, the normal component of the fluid velocity must match the normal 
component of the solid velocity (the inviscid equations allow slip in the tangential direction). Let r = r(t) denote a 
point on the surface of the body B, and n = n(r) the outward normal to the body, then 

n T r(i) = n T v(r{t),t). (15) 

In addition, the surface force per-unit-area at each point on the body is given by the local force per-unit-area exerted 
by the fluid, 

f.(r(t)) = -np(r(t),f). (16) 
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3. A partitioned FSI algorithm for the one-dimensional Euler equations and a rigid body 
added mass terms 



In the recent series of articles Q] [5] [3J, a stable interface projection scheme was developed for the problem of 
coupling a compressible fluid and a deformable elastic solid of arbitrary density. The key result from [l] [2] can 
be distilled from the consideration of a one-dimensional Riemann problem consisting of a linearized compressible 
fluid (equations |21[ ) on the right with state (po, vq, cq), and a linear elastic solid on the left with state (po, vo, &q)- 
Arguments based on characteristics were used to show that for positive times the interface values (wr, err) are given 
in terms of an impedance weighted average of the fluid and solid states, 

ZV + ZV (To — CTO /, „s 

VI = + -r— , (17) 

z + z z + z 

z~ 1 a + z~ 1 a vo Vq 

°l = — i T -i 1- -i - ( 18 ) 

z 1 + z 1 z 1 + z 1 

Here z = pc v is the solid impedance based on the speed of sound, c p , for compression waves in the solid, while z = pc 
is the fluid impedance based on the speed of sound, c, in the fluid. In [TJ [2] it was found that using a projection to 



impose (171 and (181 as interface conditions resulted in a scheme that remained stable, even in the presence of light 
solids when the traditional FSI coupling scheme fails. See [TJ [2j [3] for further details. 

The present situation of a rigid body can be considered through a limit process where c p becomes large compared 



to c, and the elastic body becomes increasingly rigid. Taking the formal limit z/z — > oo in equations (17l-(18l, with 
z fixed, results irQ 

vi = v , (19) 
oi = ao + z(v — Vo). (20) 

Thus for a rigid body, the interface surface stress is equal to the stress from the fluid plus z times the difference of the 
fluid velocity and the velocity of the body. The dependence of the interface stress, 07, on the velocity of the body, 
vo, has thus been exposed. 

These interface conditions can be derived more directly by considering the Riemann-like problem, shown in Fig. [T] 
that consists of a rigid body of mass m b adjacent to a compressible fluid governed by the linearized Euler equations. 
Using characteristic theory, we can write an explicit equation for the motion for the rigid body in terms of the initial 
conditions. This process introduces an added mass term into the equations, and the motion of the body is seen to 
be well defined even when nib = 0. The equations are then written in an alternative form as an interface projection 
that is localized in space and time. This form can be used to generalize the approach to multiple dimensions. 

Consider then the solution to the linearized one-dimensional Euler equations for an inviscid compressible fluid, 
in the moving domain x > rt(t) as shown in Fig. [T] 

dtp + vd x p + pd x v = 
d t v + vd x v - (l/p)d x a = , for x > r b (t), (21) 
d t (J + vd x (T — pcd x v = 
|>(af,0),«(ar,0),ff(x,0)] = [p (x),v (x),a (x)]. (22) 

Here a — —p is the fluid stress. The equations have been linearized about the constant state [p,v,p\. The linearized 
speed of sound is c = \J jp/p and the the initial conditions are given by [po(x), Vo(x), ao(x)]. The fluid is coupled to 
a rigid body of mass m;, whose motion is governed by Newton's law of motion for the velocity, «(,, and the position, 
Xb, of the center of mass, 

m b v b = a(n(t), t)Ab + fb, (23) 
x b = v h . (24) 

Here A b is the cross-sectional area of the body, /(, is an external body force and Tb = Xb + Wb /2 defines the point on 
the body that lies next to the fluid (wb being the constant width of the body). 

From the theory of characteristics, the variable % — a + zv i s constant along the C~ characteristic dx/dt = — s = 
v — c. Therefore, for a point Tb(t) on the body, x( r b,t) = x( r b + s^O), and thus 

a(r b ,t) + zv(r b ,t) = a (r b + st) + zv (r b + st). (25) 

Using the interface condition v(rb,t) = Vb(t) it follows that the stress on the body is 

a(n,t) = a (r b + st) + z(yo(fb + st) — v b ). (26) 



"This limit process could be quite complex and we are speaking here on informal grounds for motivational purposes. 
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C : a + zv — (Tq + zvq 



Figure 1: The x-t diagram for the one-dimensional fluid/rigid-body problem. 



Substituting ( 26 \ into ( 23 1 gives an equation for the motion of the body that only depends on the initial data in 



the fluid and the external body force, 

m b v b = a (r b + st)Ab + zA b (v (r b + st) — v b ) + f b (t), 



r b = v b . 



This equation can be written in the form, 



m b v b + zA b v b — ao{r b + st)A b + zA b va(r b + st) + f b (t), 
h = v bl 



(27) 
(28) 



(29) 
(30) 



where the added mass term zA b v b has been moved to the left-hand side. Note that equations (|29[)-(|30|l can be used 
to solve for v b even when m b = (provided zA b > 0). By using an ODE integration scheme that treats the added 
mass term zA b v b implicitly, equation ( 29 1 can be used to evolve the rigid body with a time step that need not go to 
zero as m b goes to zero. 



In practical implementation, it is often beneficial to localize (261 in space and time. Using x(i) = — e ) along 
the C~ characteristic and letting e — > leads to the relation 

a(r b , t) = a(r b +, t-) + z(v(r b +, t-) - v b (t)) . (31) 

Here a(r b + , t—) and v(r b + , t—) denote the stress and velocity in the fluid at a point which lies an infinitesimal distance 



backward along the C characteristic. Equation (31 1 is in a form that can be used in an interface projection strategy 



and can be generalized to a multidimensional problem as is done in Section [6] Furthermore, notice the similarity of 
( |31| l to equation ( |20[ ). This hints at the close connection between ( |3 1 1 ) and the projection schemes evaluated in [I] [2] [3] 
for coupling compressible fluids and deformable bodies. 



4. Normal mode stability analysis of the one-dimensional FSI model problem 



In order to understand the stability of a numerical scheme that uses the new interface conditions (31 1, consider 
the one-dimensional model problem of a rigid body confined on either side by an inviscid compressible fluid, as shown 
in Fig. [5] As in pQ we can linearize and freeze coefficients about a reference state to arrive at a problem where the 
equations of acoustics govern the two fluids, and Newtonian mechanics govern the motion of the solid. As shown in 
Fig. [2] the body has a width of w b and its cross-sectional area is assumed to be 1. Note that the equations for the 
fluids are defined in fixed reference coordinates, x < —w b /2 and x > w b /2. 

More specifically, the governing equations for the fluid in the left domain are given by 



d 


Vl 







J_- 

PL 


d 




dt 


<JL_ 




_Plc\ 


_ 


dx 





while those for the fluid in the right domain are 



d 


VR 







1 - 

PR 


d 


VR 


dt 


<JR_ 




PRCr 


_ 


dx 


or 



0, for x < 



= 0, for x > 



w b 
' 2 ' 



Wb 



The motion of the rigid body is governed by 

m b v b = T, 

where the force exerted on the rigid body by the fluid is 

T — CTR 



x = — w 



,/2 • 



(32) 

(33) 
(34) 
(35) 
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-w b /2 



w b /2 



x 
-> 



fluid (acoustics) 



solid (rigid body) 



fluid (acoustics) 



-3 -2 -1 
vl 



Vb 



vr 



Figure 2: Schematic of the one-dimensional FSI model problem used in the stability analysis. A solid rigid body is embedded 
between a fluid domain on the left and a fluid domain on the right. The boundaries of the rigid body are located mid-way 
between the ghost points of the fluid grids with index i = and the first grid point inside the domain with index i = —1 on the 
left and i = 1 on the right. 



The system is closed using interface conditions at x — ±Wb/2 which enforce continuity of velocity, namely 



VL 



Vr 



„/2 = «*» 



= Vb. 



(36) 
(37) 



Notice that the problem is posed in a moving reference frame (which we call x), and the frame attached to the rigid 
body can be calculated as 

x — x + Vb(r) dr. 
Jo 

4-.1. A first-order accurate numerical discretization of the model problem 



This section describes the discretization of the governing equations (32h(34l to first-order accuracy. As in [I], 
Godunov style upwind schemes will be used to discretize the fluid domains. We will analyze and demonstrate the 
properties of these schemes when combined with various discrete interface conditions. The finite difference grid for 
the discretization of the one-dimensional problem is outlined in Fig. [2] Note that the left and right boundaries of 
the rigid body are located at the mid-point of computational cells. This choice is made for convenience, but is not 
critical to the analysis. The grid points to the left of the rigid body are denoted by 



and to the right by 



xla = 



xr,, 



+ 



+ 



Ax L , 



Ax R , 



-2,-1,0, 



Ghost points, corresponding to index i = for both domains, will be used to enforce the interface conditions. 

Let Zk — pkCk denote the acoustic impedance in domain k = L,R. The eigen-decomposition of the matrices in 
( 32 l and ( 33 1 is given by 





Pkc'i 



Pk 




— Rk^kR 



Rk — Cfc 



-1 1 

Zk z k 



Ak = 



-Ck 
c k 



Rk L = 



2ckZ k 



-z k 1 
Zk 1 



(38) 



k,i ~ Vk(xk,i,t n ) and <rj? i w Uk(xk,i,t n ) denote discrete approximations to the velocity and stress at time 



Let u£ 

t n — nAt. We also use the notation 



Vk 


n 


v k,i 


_<?k_ 


i 


a k,i_ 



The first-order accurate upwind scheme is given by 



Vk 


n + 1 


Vk 


n 


Ok_ 


i 


Uk_ 


i 



+ AtR k A k -R^ L D- 



+ AtR k AJ Rk -D-t 



(39) 



for i = . . . , —2, —1 on the left and for i = 1, 2 ... on the right. The negative and positive parts of the wave speed 
matrices are defined by 



Ak = 



-c k 




and At = 





c fe 



respectively. The forward and backward divided difference operators are defined by D+Ui = (iti+i 
D-m — D+Ui-i, where Aa; is taken for the appropriate domain. 



(40) 

iii) /Air and 
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The methods we consider can be presented using a unified notation. Motivated by the discussion in Section [3] 
the interface stresses at t n+1 for the first-order scheme are defined by 



n + l 



n + l . / n+l n+l \ 

= <?l,-l + °i K - V L-l) 



,i +a R (v R1 



(41) 
(42) 



where oll and olr are parameters that can be used to obtain various discrete interface conditions. The traditional 
interface coupling approach found in the literature can be described in words as applying the velocity from the solid 
as a boundary condition on the fluids, and applying the stress in the fluid to derive the applied force on the body. 
This condition is achieved by setting otk = 0. Our new projection interface condition is given by setting a.% = Zh- 

The solution state in the ghost cells at t n+1 is defined to first-order accuracy by imposing continuity of the velocity 
at the interfaces 



n+l n 



and extrapolation of the stress to first-order accuracy as 



n+l 
O r n = O 



n+l 

a B n = a 



n+l 
I,L i 

n + l 
I,R ■ 



The rigid body equations (341 are advanced in time with the backward Euler scheme, 



m b v b 



rribv b + AtJ 7 



where the force at t n+ , J rn+1 , is defined as 



J 7 " 



(43) 
(44) 

(45) 
(46) 

(47) 

(48) 



The backward Euler method is used here in order to simplify the analysis. Used in isolation, the backward-Euler 
scheme is unconditionally stable for any At independent of m b provided m b > 0. We will show, however, that 
the fully coupled FSI problem has a time-step restriction that depends on m b for the traditional interface coupling 
scheme. For the new interface projection scheme we show that there is no dependence of the stable time step on m b . 
The backward-Euler scheme is, of course, only first-order accurate. For higher-order accuracy one can use implicit 
Runge-Kutta schemes, as described in Section [7] where we extend the scheme to multiple space dimensions. Note 
that while implicit schemes may be more expensive per time-step than explicit schemes, they are only used to solve 
the rigid body equations which consist of just a few ODEs. As an alternative to implicit schemes, one can consider 
using an explicit scheme with a sub-cycling approach (i.e. taking multiple sub-steps with a smaller value for At). 
Some remarks on these issues will be provided in subsequent discussions. 

In summary, to advance one time level from t n to t n+1 using the first-order accurate scheme, the following steps 
can be followed 

Algorithm 1. 



1. Compute 



1 n+l 



for i 



1 and 



Vr 

OR 



-i n+l 



fori = 1,2,... by 



39) 



2. Set T n+L = tr£7i + QflK,i ~ < ) - ~ « ~ <t_i), and solve (47) for 



„n + l 



v£ +1 = \m b + At(a L + a R )| \m b v b + At^cr^ 1 + a R v^ - (o^t^i - aL«2 + _\)j] . (49) 
3. Define the ghost point values by the velocity boundary conditions \43ty and \44\> along with the stress extrapo- 



lations (45) and U 



4-2. Normal mode analysis of the first-order accurate scheme 

Next, we analyze the stability of the interface discretizations, and investigate how the choice of a_t and an affect 
the behavior of the overall numerical method. To simplify the presentation, assume cl = cr = c, pl — Pr — p, 
Axl ~ Axr = Aa;, and a_L = ctR = a. In addition set z — zl — zr. These assumptions are purely for convenience 
and clarity, and do not materially change the results of the analysis. We pursue a stability analysis via the normal 
mode ansatz of Gustafsson Kreiss and Sundstrom [B]. 

As was done in [l], we seek normal mode solutions of the form 

v b = A n v b , for k = L, R, (50) 
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where Vk,i and <Jk,i are bounded functions of space, and A, the amplification factor, is a complex scalar with \A\ > f . 
If such a non-zero solution can be found (for given values of the parameters A, At, mj, z, etc.) then there are solutions 
that grow in time and we say that the scheme is unstable for those parameters. We note that more general definitions 
of stability allow some bounded growth in time, but for our purposes here we use this more restrictive definition. 
Characteristic normal modes are denoted by 



Rl 



1 

2~cz 



o k — zv k 
CTjfe + zv k 



(51) 



for i 



-3,-2,-1 



for 4 = 1,2, 3, 



(52) 



(53) 



Insertion of ( 50 1 into the finite difference scheme ( 39 \ leads to 

AaL,i = at, t i — A {ah,% — ac,i-i) 
Mh,i = b L ,i + XQ>L,i+i - bt,i) 

and 

Aclrj — an,i — A (oij.j — or^-i 

AbR t i = bn,i + A (6fl,i+i — 6n,i) 
where < A = cAt/Ax < 1. Define the quantity 

A-l + X 

We see that \r\ > 1 by rewriting jrj 2 > 1 in terms of the polar variables R and 6, where A = Re 10 . By simple 
algebraic manipulations, |r| 2 > 1 can be rewritten as 

{R - l) 2 + 2\(R - 1) + 272(1 - A)(l - cos0) > 0, 

which is true since R > 1 and A < 1. 



For the two components on characteristics coming in from infinity, the solution to the difference equations ( 52 1 

-2,-1, 



and ( 53 \ is 



ah, 



r 'a.L,-i, 



for i — . . . , —3, 
for i = 1,2,3,.. 



The assumption of boundedness as i — > ±oo gives a,L,i = for i . . . , —3, —2, —1, and bn,i — for i 
that cll,o and &ij,o do not play a role in the difference equations ( |52[ ) and (J53b, but their values can be determined 
algebraically using the interface conditions 



a-L.o 
bR,o 



2z 



2z 



-{vb/c— b L ,~i) 
(v fl /c + a Ril ). 



The remainder of the solution to difference equations ( 52 1 and ( 53 1 is given by 



aR,i 
bh.i 



r Or,o, 



for i 
for i 



0,1,2,3,..., 
...,-3,-2, 



-1,0 



(54) 
(55) 



The solutions (54 1 and (55 I are bounded because jr| > 1. The definition of the characteristic normal modes on the 
interior yields 

[%-,! Til . 

r l bL,o for i = . . . , —3, - 



V 




T 




= c 




a 




z 




L,i 





-1 



and 



V 










— c 




r 


a 


R.,i 













for i = l,2, 3, 



(56) 
(57) 



The three undetermined constants 6l,o, clr.o, and Vb are defined by application of the interface conditions ( 43 1-( 46 1 
and the rigid body integrator (47 1. This leads to the linear system of equations 

z 



1 + 



2zr 







AAt 



(z - a) 



1 + - 
AAt 



2zr 



(z~a) m b (A ~ 1) + 2AtaA 



a + z 






2z 




bh.a 


a + z 






2z 




Vb/c_ 



(58) 



The system (581 is an eigenvalue problem for A, in the sense that if there is an A such that the determinant of the 
system is zero, then there exists a non-trivial solution of the form If, furthermore, |^4| > 1, then the solution 

(50 1 grows in time. 
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Theorem 1. The numerical scheme using the interior discretizations \39\j , interface conditions {4$j -\ rigid body 
integrator (41) and projections 41 )-\4%) with a = z has no eigenvalues A with \A\ > 1 for A < 1 and m;, > 0. 



Proof. For a = z the eigenvalue problem 1 58 1 reduces to 



= 0. 



a-R,o 

Vb/c_ 

The determinant is zero when A = mb/(mb + 2Atz). By assumption, At > and z > and so \A\ < 1. 



1 -1 

1 1 

rn b (A - I) + 2AtzA 



(59) 



□ 



Theorem 2. The numerical scheme using the interior discretizations (39), interface conditions \43\ -\4d), rigid body 
integrator (^7| ) and projections \4 l\ -\4^ with a = has no eigenvalues A with \A\ > 1 when 

At < m 6 (4- X)/(z\) (60) 

for A < 1. Conversely, if At > m{,(4 — A)/ (zX), then there are eigenvalues with \A\ > 1 for A < 1. 



Proof. For a — 0, the eigenvalue problem (58 1 reduces to 

1 



1 2r 


AzAt 





_4zAt 



1 

1 

2 



m 6 (A- 1) 



The zero determinant condition is solved to give three roots Ai = 1 — A/2 and 



A 2 ,3 = 1 



A z£A 



4 2 

where £ = At /mi,. Clearly, |^4i| < 1 for A < 1. In the case 



A z£X 
4 2 



A z£\ 



1 + 



A2 and ^3 are complex conjugate, and 



|^ 2 | 2 = |^3| 2 = (1- J- 



A z£A 

2~ 



2 _ fi _ A _ ^ 2 
4 2 



+ 1 = 1 < 1. 

2 2 



When A2 and ^3 are real, rewriting ( 62 1 as 



. .A z£X. , 

A 2 ,s = l-(^ + -|-)± 



A z£A 
4 + ^T 



shows directly that both roots are are always < 1, hence |«4| < 1 if and only if 



A z£X 
4+^ 



which is equivalent to 



A rfA 
4 2 



The necessary condition that the right hand side is positive is equivalent to 

Assume ( |64[ ) holds and square both sides of ( |63[ ) to obtain 

z£A < 4 - A. 



(61) 



(62) 



Hence, \A\ < 1 exactly when (651 holds. The proof is completed by observing that (65 I is equivalent to (60 1. 



(63) 
(64) 

(65) 

□ 
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Remark: Theorem [2] indicates the traditional coupling scheme with a = has a time step restriction that can 
be more strict than that for the fluid domains alone. Even though the rigid body is formally integrated with the 
Backward Euler scheme, which would be unconditionally stable when used in isolation, the coupled formulation does 
not include the full dependence of the forcing J rn+1 on v^ +1 . In the case of light bodies, i.e., bodies with small mj, 
the time step restriction for the scheme with a — 0, can be severe. Another way to state the result is that for any 
fixed grid resolution, there exists some sufficiently small mass for which the solution will have exponential growth in 
time. In fact, it is easy to see that for a = and fixed At, the limit of small mass yields r ~ 1 — A/4 — zXAt/rrib and 
therefore lim mb ^o \A\ = oo. Note, however, that the traditional coupling scheme with a = is stable, for any finite 
mass nib, provided the time step satisfies the conditions given in 2\ 

Remark: Theorem [I] shows that the time step restriction (60| can be avoided by switching to the interface 
conditions with a = z. 

Remark: The structure of the eigenvalue problem in the proof of Theorem [I] suggests why the choice a — z is in 
some sense optimal. When a = Z, the rigid body mode is decoupled from the fluid modes and stability follows for 
any mt > 0. On the other hand, for a — 0, the eigenvalue problem (58 1 represents a coupled system and the question 
of stability is summarized in Theorem [2] 

Remark: For choices of a other than zero or z, the stability of the numerical scheme varies somewhat. The 
determinant condition can be used to produce an expression for A, but it is somewhat difficult to interpret. We 
provide no further discussion about other choices of the parameter a. 

4--S. A second-order accurate numerical discretization of the model problem 

We look now at the formulation and stability of a second-order accurate version of the projection interface scheme. 
For the discretization of the fluid domains we choose the second-order accurate Lax-Wendroff scheme, 



Vk 


n + 1 


Vk 


n 


Ok_ 


i 




i 



+ AtC k D 



+ 



-ClD+D. 



k — L,R, 



(66) 



where Do = (-D+ + D )/2 is the centered difference operator, and Ck has been defined previously in ( 38 1 . The 
Lax-Wendroff scheme is a good model, since many non-linear schemes of TVD type are designed to approximate 
the Lax-Wendroff scheme in the parts of the computational domain where the solution is smooth. The projection 
coupling conditions can be implemented to second-order accuracy as follows. Define interface stresses on the left and 
right at any time t" by 



<?I,R 



+ OLL \V b 



3«2,-i - 1>£,-2 



J(7 R.l — a R,2 , I JV J{,1 — V R,2 

- + OCR 



(67) 
(68) 



2 V 2 

These are obtained by extrapolation from domain interiors, and subsequent projection. The force at any time level 



t n is defined as before using (481, and a second-order accurate trapezoidal integration for the solid is then defined 

1 



m b - 



At 



(69) 



The velocity from the solid is applied as a boundary condition on the fluids to second-order accuracy by setting the 
average (vl,o + «l,-i)/2 equal to v b (and similarly at the right interface), or equivalently 



n r» n 

V L ,0 = 2 «6 
n n n 

v R ,o = 2«i, 



Vl,-i, 

71 

Vr,i- 



Extrapolation of the stress to the ghost cells gives 



n n Tl n 

a L,0 = ^I,L - &L,-1, 

°\R,0 — Za I,R ~ a R,l- 



(70) 
(71) 



(72) 
(73) 



In summary, to advance one time level from t n to t n+1 using the second-order accurate scheme, the following 
steps can be followed 

Algorithm 2. 



1. Compute 



1 n+l 



for i 



-1 and 



vr 
or 



-i n+l 



for z = l,2,. 



using 



66) 
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2. Define F n using the computed solution at t and solve i69j) to yield 



m b - 



AtctL AtctR 



Ata L 



Ata r 



v b + 



. / o 11 + 1 

At / 3cr R1 - a 



/V 2 ■J (T if.,l — a R,2 



a R At ( 3w£+ 



2 

3<j 



At / 3^+\ 



n+l „ n _n 

'L,-2 •J <t L,-1 — "L,-2 



n+l 

Z.,-1 



2 

3ttf 



+ 



(74) 



3. Define the ghost point values for the fluid domains using (70) - \7j\ 



4-4- Normal mode analysis of the second-order accurate scheme 

For the analysis, we make the same assumptions as in Sec. |4.2[ that the grid spacings and wave speeds are the 



same on both sides of the body and that ctL = olr = a. First, we decompose (661 into characteristic components, 
and obtain the two scalar equations on each side of the body, 



c 2 At 2 



n+l n 

a k,i — a fe,i 

where k — L, R, i — . . 

a™ = A n r* and 6" = A n r l into (751. This leads to the characteristic equation 

\(v + v 2 )r 2 + (1 - A - v 2 )r + \{v 2 



+ cAtD bti 



-D+D-b" k 



(75) 



for k = R. The normal modes are found by inserting 
v) = 0, (76) 



where v = cAtj Ax for the 6 characteristic component and v = — cAt/ Ax for the a characteristic component. The 
assumption c = cl = cr gives the same characteristic equation on either side of the body. There are four roots, two 
for the — c characteristic, that we denote rf/ and and two roots for the c characteristic, that we denote rf and 
rj. It is well-known, see e.g., [6], that for the equation ut = cu x under the CFL-condition A < 1, the two roots of 
1761 satisfy 



for some 8 > when c > 0, and 



\rt\<l-S 
\rt\ > 1 
rt = 1, 



|rr| < 1 
r\ = 1, 
|ra | > 1 + 5 



\A\ > 1, 

1-4.) > 1, A^l, 
A = l, 



\A\>1, A^l, 
A = l, 



(77) 



(78) 



\A\ > 1, 



when c < 0. For the model problem ( 66 1 , there are thus four roots. From ( 77 1 , ( 78 1 and the condition of boundedness 
at infinity, it follows that the rf and r-f components are zero for i < and that the and r^ components are zero 
for i > 0. Hence, the normal mode solutions to the left and to the right of the body can be written 



and 



respectively. 



D 








= c 




a 




z 







(r 2 Yatfi + i 



(n ) l a B , + c 



(rtYb R ,o 



for i < 0, 



for i > 0, 



(79) 



(80) 



The solutions (791 and (80 1 inserted into the interface conditions (70 1, (71 1, (72 l, and (73 1 together with 1 74 1 
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give five equations for the five unknowns a R} o, b Rl o, o,l,o, 6l,o, and Vb- Fully written out these equations are 



1 + 



a-L,o 



1 + I 6i,o + -Vb = 0, 



(81) 



(l + r 1 ) a R ,o - (l + rf ) b R>0 + -«6 = 0, 

(82) 



-2 - + 1 



3J_ 

2rZ 



1 1 



1 + Ti 



2 2 (,,-)> 



H t ,o + 1+ -r +2 1 



3 1 I I 

T 

2 



2 r+ 2 (r+) 2 



&L,0 



2a t);, 



+ 1 + rf - 2 - + 1 



+ 



2a Vb 



0. 

(83) 
0, 

(84) 



3 1 

2^ 



1 1 

2(^F 



a\/3 1 11 



2 r + 2 ( r +)2 



+ 



&L,0 



A — 1 luib 2a\ Vb 
A+l Atz + z ) c 



For the case a — z the system (81 1-(85 I becomes 



1 + 



a.L,a 



1 + -4- I 6i,o + -«(, = 0, 



v l + r x ) a fl>0 - (1+ r{) b R , + -v b = 0, 



+ 



(r 2 - 



&£,0 «6 = 0, 

C 



>1 + r-n a fl , + (1 - 5r+ + 2(r+) 2 ) 6h,o + -«(, = 



( 3 1 "\ Cq + / +^2^ i. , (A - 1 2m b \ v b 

Irf - (r^J ^° ~ (3ri " (ri ) } **° + [a+1 Mz + 2 ) c = °- 



(85) 

(86) 
(87) 
(88) 
(89) 
(90) 



Theorem 3. When \A\ > 1 and A < 1, t/ie system '86)- '90) only has the trivial solution o,l,o = 6l,o = or,o = 
°r,o — v b = 0. Hence, the numerical scheme using the interior discretizations \6(fy, interface conditions (70|)-(73|), 



rigid body integrator (74-) and projections 61)-\68^ has no exponentially growing modes for A < 1 and m,b > 0. 



Proof. Adding equations ( 86 1 and ( 88 1 gives 



2 1 



ol,o = 0. 



Because of ( 78 \ , 



r~ ~ \r7 ~ 1 + 5 1 + 5 



>0, 



'2 l'2 

and consequently, a^.o = 0. Similarly, subtracting ( |87[ ) from ( |89[ ) and using ( |77[ ) gives fo^o = 0. Equation ( |90[ ) with 
ai,o = 6jj,o = gives 



A non-trivial solution exists if 



A - 1 2m 6 
.4+1 Ate 

.4-1 2m b 



+ . ) o. 



2 = 0, 



(91) 



.4+1 Atz 

which is equivalent to A = (mt — Atz)/(m,b + Atz). Assuming that for At > 0, z > 0, and mi, > 0, it follows that 
\A\ < 1, and hence that the only solution of (91 1 when |.4| > 1 is Vb = 0. Finally, the remaining equations 

1 



b L .o = 



and 



(1 + r x )or, = 0, 



have the unique solutions 6l,o = or,o ~ 0, because (771 and (781 exclude the possibility that = — 1 or r x = — 1 
when L4I > 1. □ 
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t Ax 

Figure 3: Results for the one-dimensional FSI problem with mj, = 1 for the first- and second-order accurate schemes. Top 
left: velocity at t = 0.75. Top-right: stress at t = 0.75. Bottom left: velocity of the rigid body, vt, versus time. Bottom right: 
convergence of the max- norms errors (reference lines of the corresponding order are displayed in black). The solutions are 
plotted in the reference domain [— f , 0] for the left domain and [0, 1] for the right domain with the rigid body of width uij, = 
located at x = 0. 



5. Numerical demonstration of the theory for the FSI model problem 

We now present numerical results from solving the one-dimensional FSI problem introduced in Section [3] The 
aim is to demonstrate the accuracy and stability of the new FSI projection algorithm For this purpose we use the 
exact solution derived in | Appendix A| The problem consists of an initial Gaussian pulse in the fluid that moves left 
to right and interacts with the rigid body. The initial conditions for the velocity and stress are given by 

v(x,t = 0) = ~ exp (~/3 2 (x - x ) 2 ) , a(x,t = 0) = ex P {~P 2 { X ~ x o) 2 ) ■ 



The rigid body is initially at rest. The exact solution is denned by (A. 191, (A.7l and ( |A.10 1. Throughout this section 



we use pl = 1, cl = v2, Pr = 1, cr = v3, /3 = 10 and x<j = —1/2. Note that the initial conditions (A. 171 and 
( |A.18[ ), and exact solutions ( |A.19| |, ( |A.7[ ) and ( |A.10[ ) may require differentiation with respect to space and/or time 
in order to be used or compared with the dependent variables of velocity and stress which we use. 

5.1. Easy case: rigid body with mass one 

We begin our numerical results with a case where the CFL time-step constraint in the fluids is dominant over the 
explicit ODE time-step constraint for the rigid body. This is the case when 



max(zi,zji) < m^mm 



( Cl Cr \ 
\Ax L ' Ax R I ' 



which implies that time steps which satisfy the usual CFL stability constraint in the fluid also satisfy the stability 
constraint associated with the ODE for rigid body motion. As a result, the traditional interface coupling technique 
found in the literature has no difficulty, and we are simply setting out to demonstrate that the new interface projection 
technique remains accurate for this case. 
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Figure [3] shows simulation results for m& = 1 when using the first-order accurate upwind scheme for the two 
fluid domains, the backward Euler integrator for the rigid body evolution equation, and the interface projection 
scheme with a = z as defined by Algorithm [T] In addition we show results using the second-order accurate Lax- 
Wendroff scheme for the fluid domains together with the trapezoidal rule for integration of the rigid body as defined 
by Algorithm [2] with a — z. 

For both cases we use Axl = Axr = 1/50. The exact solution and numerical approximations for v and a are 
displayed as functions of the reference coordinate x at t = 0.75, and the velocity of the rigid body is shown as a 
function of time. The width of the body is taken as % = (this has no influence on the results) so that the left and 
right reference domains meet at x — 0. The results from the first-order accurate scheme show predictably smeared 
out solution profiles. The results from the second-order accurate scheme are in very good agreement with the exact 
solution even at this coarse resolution. Figure [3] also presents results from a grid convergence study and shows 
the max-norm errors for this problem using the two algorithms. The predicted convergence rates are convincingly 
demonstrated for both velocity and stress. 

Remark: For this case, one could also use the new projection scheme with a forward Euler rigid body integrator. 
Simulation results for this case reveal no unexpected behavior. 

Remark: For this case, traditional coupling techniques without projection would not experience exponential 
blowup for the considered grids and time steps. Numerical results using the traditional scheme with a — for this 
case are nearly identical to those in Figure [3] and are therefore not shown. 

5.2. Difficult case: very light rigid body with mass 10~ 6 

We now consider a case where the time-step restriction for the traditional interface algorithm is orders of magnitude 
smaller than the time-step restriction for the new interface projection algorithm. The time-step restriction for the 
new projection algorithm depends only on the usual CFL time-step restrictions for each fluid domain separately; the 
coupling with the rigid body imposes no new constraint on the time-step since the backward Euler and trapezoidal 
methods are both A-stable. We consider a rigid body with mass mi = 10 -6 and use the same grid spacings as before, 
Axl = Axr = 1/50. Figure |4] shows simulation results for this case using the two new schemes. As the figure shows, 
the results from the second-order accurate scheme are, as expected, superior to those from the first-order accurate 
scheme. The lower right graph in Figure [4] presents a convergence study. The expected rates of convergence are again 
convincingly demonstrated. 

Remark: For this case, one could instead consider using an explicit rigid body integrator together with the 
projection scheme. The rigid body integration must respect the ODE time-step constraint and so subcycling can be 
used. It is straightforward to estimate that for Axl = Axr — 1/50, 14389 subcycles are required to obtain stability 
of a forward Euler integrator. The number of subcycles required for stability decreases, however, as At decreases. As 
a result, for Axl = Axr = 1/1280 (the finest resolution in the associated convergence studies), only 568 subcycles are 
required. A sub-cycling has been implemented and the results are nearly identical to the results shown in Figure [4] 

Remark: Had the traditional algorithm with a — been used, the entire solution (both fluid domains and the 



solid domain) would have to be integrated using a time-step which satisfies a constraint of the form of ( 60 1 . For the 



first-order scheme with Backward Euler rigid body integrator the constraint is (60 1. For other fluid discretizations 
and/or rigid body integrators,the timestep restriction can be determined following the approach used in the proof of 
Theorem [2] Such a time-step restriction can be quite severe and arises as a result of using a partitioned algorithm 
without the interface projection. 

5.3. Rigid body with zero mass 

The new projection based FSI scheme remains well defined even when the mass of the rigid body is zero. This is 
apparent from the update equation for the velocity of the rigid body, equation Q49J) for the first-order accurate scheme 



or equation ( 74 1 for the second-order accurate scheme. The traditional partitioned algorithm is not well-defined for 
this case, since it would require division by m&, and so is not an option. Figure [5] shows results for the first-order 
upwind method with backward Euler rigid body integration, and the second-order upwind method with trapezoidal 
rigid body integration. The exact solution is computed for mi, = which yields essentially the same solution used for 
the two domain model problem in pQ, i.e. the solution behaves as if the rigid body were not present. Figure [5] shows 
convergence results where again the pred icted rates of convergence are demonstrated. No significant differences from 
the mass m& = 10 -6 case in Section [5^2] are observed. 

Remark: For this case it is impossible to satisfy the ODE stability constraint without using an A-stable integrator 
and so explicit rigid body integration with subcycling is not an option. Put another way, the explicit algorithm would 
require an infinite number of subcycles. 



6. The multi-dimensional interface approximation and added-mass matrices 

In this section we extend the added-mass algorithm to multiple space dimensions. Formula ( |31[ ) relates the 
pressure and velocity of a point on the body to the nearby pressure and velocity in the fluid. This relation is in a 
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t Ax 

Figure 4: Results for the one-dimensional FSI problem with m b = 10 -6 for the first- and second-order accurate schemes. Top 
left: velocity at t = 0.75. Top-right: stress at t = 0.75. Bottom left: velocity of the rigid body, v b versus time. Bottom right: 
convergence of the max-norms errors. The solutions are plotted in the reference domain [—1,0] for the left domain and [0, 1] 
for the right domain with the rigid body of width w b = located at x = 0. 



form amenable to multidimensional generalization. Let r = r(t) denote a point on the surface of the body B, and 
n = n(r, t) the outward normal to the body, then in multiple space dimensions (311 becomes 

-p(r(t), t) n = -p(r+, t-)n + z(r+, t-) [n T (v(r+, t-) - v(r, t))] n, 

where v(r, t) — r is the velocity of the point. To clarify the notation let p r = p(r, t) and v r = v(r, t) denote the 
pressure and velocity on the body at point r = r(t), and zj = z(r+, t—), Pf = p(r+, t—) and v/ = v(r+, t—) denote 
the impedance, pressure and velocity at the adjacent point in the fluid. This gives 

— p r n = — p/n + zf [n T (v/ - v r )] n. 



Using equation |13| for v r = r it follows that 

— p r n = — pfn + zf [n T (v/ — V6 + Yu>)] n. 



The key point of ( 92 l is that is shows how the force exerted by the fluid on the body, f s 



(92) 

-pr-n, depends on the 



velocity of the center of mass, Vf,, and the angular velocity, lj, of the body. Substituting ( 92 1 into the expressions ( 10 1- 
( 11 1 for F and T gives 

T = / z/nn T (— v b + Ylj) ds + / -p/n + 2:/(n T v/)n ds + f b , 

JdB JdB 

T= / ZfYnn T (-v b + Ylj) ds + / y x ( — p/n + z; (n T v/)n) ds + g b . 

JdB JdB 



We write T and T in the form 



F = -A m v b - A vu u + J=, 



T 



A OJV AbJU) . /-I— 

-A v b -A uj + T 
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Figure 5: Results for the one-dimensional FSI problem with mj = for the first- and second-order accurate schemes. Top 
left: velocity at t = 0.75. Top-right: stress at t = 0.75. Bottom left: velocity of the rigid body, versus time. Bottom right: 
convergence of the max-norms errors. The solutions are plotted in the reference domain [—1,0] for the left domain and [0, 1] 
for the right domain with the rigid body of width uij, = located at x = 0. 



where the added-mass matrices A 1 ^ are given by (using Y T — —Y, where Y is denned by ( 14 1 ) , 



2/nn T ds, 



A v 



f z f (Yn 
J as 



J dB 



ds, 



n 1 ds A" u = I z f Yn(Yny ds, 

dB 



and T and T are given by 



-p/n + z/(n v/)nds + ff, 



T= I yx (- p f n + Zf(n v/)n) ds + g b . 

Jbb 



(93) 
(94) 

(95) 
(96) 



Note that A vv and A"" are symmetric and positive semi-definite while (A VUJ ) T = A" v . Let A m € K 6x6 denote the 
composite added mass matrix (tensor), 



Am, — 



A vv A vu ~ 
A uv A u " 



This matrix is symmetric and positive semi-definite since for any vector w = [a b] T G R 6 , a £ R 3 , b G 



w A m w = [a b] 



A vv A vul 

j^v A uiui 



z f I ||n T af + 2(n T a)(yn) T b + \\(Yn) T h\\ 2 )ds, 



2 / ((n T a) + (yn) T b) ds. 
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The rigid body equations of motion |5|-([8]l can now be written in the form 
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(97) 



We will refer to equations (97 1 as the added-mass Newton-Euler equations. 

Remark: By solving equations (971 with an implicit time-stepping scheme that treats the added-mass terms 



implicitly, the rigid body equations can be advanced with a large time step even as m;, and A approach zero, provided 
Am is nonsingular. This is described in more detail in Section [7] 

Remark: In | Appendix B| we present the form of the added-mass matrices for some common body shapes. 



7. The multi-dimensional time-stepping algorithm 

We make use of overlapping grids to treat multi-dimensional problems with moving rigid bodies. Narrow boundary 
fitted grids lie next to the bodies and these move with the bodies (see the examples in Section [8| . One or more 
stationary background grids generally cover most of the domain. This approach results in high-quality grids even as 
bodies undergo large motions. The time-stepping algorithm we use for FSI problems with rigid bodies is described in 
detail in [7J, while that for FSI problems with deforming solids is described in [2]. In [7] the Newton-Euler equations 
for the rigid bodies are solved using a Leap-frog predictor step followed by a trapezoidal rule corrector step. 



The FSI time stepping algorithm 
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Figure 6: The FSI time stepping algorithm for advancing the states of the fluid and rigid body. 



For the new interface algorithm developed here, we choose a time-stepping method for the Newton-Euler equa- 
tions (971 that treats the added-mass terms implicitly so that the scheme is well-defined even in the limit of zero 
mass and/or moments of inertia. We use diagonally implicit Runge-Kutta (DIRK) schemes for this purpose 20 . 
DIRK schemes have very nice stability and accuracy properties. The one-stage, first-order accurate DIRK scheme, 
which we denote by DIRK1, is just the backward-Euler scheme. For the numerical results in section [8] we will use 
a two-stage third-order accurate (A-stable) scheme, denoted by DIRK3, due to Crouzeiux (see [2D] (2.2)). In each 
stage of the DIRK scheme we solve an implicit approximation to (971 by Newton's method. 




Figure 7: The fluid grids for two-dimensional problems have a grid point aligned with the boundary of the rigid body. The 
solution on the boundary is wj, while w™ 2 an d denote the values on the ghost points. For clarity, only one grid line is 

shown in the direction normal to the boundary. 

The FSI time stepping algorithm for advancing the fluid and rigid body is outlined in Figure|6] In a slight difference 
from the grid arrangement used for the analysis in one-dimension as illustrated in Fig[2] the two-dimensional grids 
have a grid point aligned with the boundary of the body as shown in Fig[7| Let w" = (p™, v",p") denote the discrete 
solution in space and time for the state of the fluid at time t", where i is a multi-index. Let (x" , v^ 1 , oj", E n ) denote 
the discrete approximation in time to the state of the rigid body. Let x™ = denote the (moving) grid points on 
the fluid grid that lies next to the body and G" the grid velocity (the fluid domain will actually be discretized with 
multiple overlapping grids but for clarity we ignore these other grids in the current discussion) . 
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Suppose that we are given the full state of the discrete solution at time i" -1 and wish to determine the state 
at the next time step t n . In the first stage of the time stepping algorithm, predicted values are obtained for the 
state of the solid body at the new time, (x£ , v p , u> p , E p ). These values can be obtained either from the Newton-Euler 
equations of motion or using extrapolation in time (for a second order accurate scheme we extrapolate using the 
current level and two previous time levels). From the predicted state of the body we can obtain predicted values for 
the grid location, G?, and grid velocity, G?; these values are needed to advance the fluid state. Note that the grids 
move as a rigid body and do not deform. In the second stage of the time stepping algorithm we obtain the new values 
of the fluid state wj 1 = (p™, v",^) at interior grid points, i G X/, and predicted values, w? = (p? , v?,p?), at points 
on the body surface, i G Tb- These values are obtained from our high-order Godunov-based upwind scheme [7]. Note 
that at this stage we do not impose any boundary conditions on the fluid at the body surface. Predicted values of the 
fluid on the body, i G Tb, can be obtained, however, since the fluid state contains ghost points, i G Xg (two lines of 
ghost points being used for our second-order Godunov scheme). Given the predicted fluid states w? we can compute 



the partial body forces ( 95 l-( 96 1 and the added mass matrices (93l-(93l using numerical integration over the surface 
of the body. We then solve the added-mass Newton-Euler equations (971 (e.g. with a DIRK scheme) to determine 
the corrected state of the rigid-body at the new time, (x" , , u> n , E"). The predicted state of the fluid on the solid 
body is then corrected by setting the fluid velocity equal to the (local) body velocity and the fluid pressure to equal 
its projected value, 

vr = v£i, iGXs, (98) 

-pr = -P? + * p n T (vf-v?, 1 ), iGX s . (99) 

Here the local body velocity is V(,,i = v^ 1 + W n (r" — x"; ), where r™ denotes the location of a point on the body surface, 
and where W n is defined from <jj n using Q. After projecting the pressure, the density is corrected using 

ft =P\( P -/P\) lh , iGXs, (100) 

which ensures that the entropy of the predicted state equals that of the corrected state. The fluid grid, Gp, and grid 
velocity, G", at the new time are corrected from the predicted values to match the new state of the rigid body. In 
the final stage of the time stepping algorithm, the ghost values of fluid state that lie adjacent to the body surface are 
updated using the appropriate boundary conditions and compatibility conditions, see [7J [5] for more details. 



8. Numerical results in two space dimensions 

In this section we present numerical results in two-dimensions that demonstrate the accuracy and stability of the 
added-mass interface algorithm when applied to light rigid bodies. A pressure driven light piston problem is used to 
examine the accuracy of the two-dimensional added-mass algorithm for an FSI problem with an analytic solution. A 
smoothly accelerated light rigid body in the shape of an ellipse is used to evaluate the scheme for a two-dimensional 
problem that includes the rotational added-mass terms. Solutions using the new added-mass algorithm are compared 
to the old algorithm, which is necessarily run at a small CFL number to avoid exponential blowup. Although the 
exact solution to this problem is not known, a posteriori estimates of the errors are determined from solutions on 
a sequence of grids of increasing resolution. As a final example we simulate a Mach 2 shock impacting an ellipse 
of zero mass and zero moment of inertia. This case demonstrates the robustness of the added-mass algorithm for a 
very difficult situation. Solutions to this shock driven ellipse problem are computed at varying grid resolutions and 
compared. These results include computations that use dynamic adaptive mesh refinement (AMR). 

8.1. Pressure driven light piston 

The geometry of the one-dimensional pressure driven piston problem is shown in Fig. [8] A compressible fluid 
occupying the region x > G(t) lies adjacent to a piston of mass mt and cross-sectional area Ab- The face of the 
piston that lies next to the fluid follows the curve x = G(t) as time evolves. A body force fb(t) also acts on the 
piston. The exact solution to this problem can be determined for a fluid that is initially at rest and the form of 
this solution is given in [7]. When fb(t) = 0, the exact solution can be determined explicitly. For general fb{t), the 
case considered here, the exact solution can be accurately approximated by numerical integration of the appropriate 
ordinary differential equations. 

We solve the pressure driven piston problem on a two-dimensional overlapping grid denoted by Qp \ where j 
denotes the grid resolution (see Figure[8|. The grid spacing in the x-direction is chosen to be Ax^ = l/(10j). The 
spacing in the y-direction is held fixed at Ay = 2/10. A background Cartesian grid covers the domain [—0.5, 1.5] x [0, 1] 
and remains stationary. A second Cartesian grid initially covers the domain [0,0.5] x [0,1] and moves over time 
according the piston motion. 

The pressure driven piston problem is solved for a piston of mass mb = 10 -6 . The initial conditions for the fluid 
are (po,«o,f>o) = (1-4,0., 1) with 7 = 1.4. The body force is chosen to be fb(t) = poAb{^ — fi 3 ) which results in a 
piston that smoothly recedes to the left and for which we expect the numerical solution to be second-order accurate 
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x = G(t) 




^piston faccjj 
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1.5 



(2) 

Figure 8: Left: the x-t diagram for the pressure driven piston problem with a receding piston. Right: overlapping grid Q p for 
the fluid region at t = 0.0. The green grid moves with the piston. The blue background grid does not move. The interpolation 
points are marked as black dots. 



Pressure driven light piston, M=10 
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Figure 9: Results for a pressure driven light piston of mass mj = 10 . Left: computed and exact solution at t = 1. using Q p 
Right: maximum errors and estimated convergence rates at time t = 1. 



(8) 



in the max-norm. The computed and exact solutions are shown in Fig. [9] for results using grid Q p 8 ^ and these are 
in excellent agreement. Figure [9] also gives the max-norm errors for solutions computed on a sequence of grids of 
increasing resolution. The values in the columns labelled "r" give the ratio of the error on the current grid to that on 
the previous coarser grid, a ratio of 4 being expected for a second-order accurate method. The convergence rate, /3, 
is estimated from a least-squares fit to the log of the error equation e(h) = ChP . The results show that the solution 
is converging at close to second-order. 



8.2. Smoothly accelerated ellipse 

In this example we consider a light rigid body in the shape of an ellipse that is accelerated by a smoothly varying 
body force. We compare the solution from the new added-mass algorithm to that from the old algorithm, the latter 
requiring a very small time step to avoid exponential blowup when the mass of the body is small. 



The overlapping grid for this rotated- ellipse problem is denoted by Qrl where j denotes the grid resolution (grid 



Qre^ is shown in Fig. 



|10l 

with grid spacing As u> — l/(10j) 



The grid consists on a stationary background Cartesian grid for the region [—2, 2] x [—2, 2], 
A narrow boundary fitted grid is located next to the surface of the elliptical body, 
and this grid will move to follow the motion of the body. The surface of the body is defined by an ellipse, which 
has major and minor axes of lengths 1.4 and 0.7, respectively, and which is rotated by 7r/4 in the counterclockwise 
direction. The boundary fitted grid extends 8 grid lines in the normal direction (the grid in Fig. 10 shows an additional 



ghost line), and the grid spacing in the normal direction is slightly clustered near the ellipse surface. The number of 
points in the tangential direction is chosen so the grid spacing is approximately As"'. 

The ellipse is accelerated using a body force that smoothly ramps from zero to one on the time interval [0, |] and 
then smoothly ramps back to zero over the interval [|, 1], In particular, the body force is in the ^-direction and is 
given by 



U(t) = R{2t) - R(2t - 1), where, R(t) = < 





35t 4 
1 



84r + 7or 



20f 



if t < 

if < t < 1 • 

if t > 



(101) 



20 



Accelerated ellipse: M=0.001 
0.8, , , , 




Figure 10: Accelerated ellipse. Left: overlapping grid C?re' at time t = 0. Right: time histories of the rigid body velocity (vi, 1*2)1 
angular momentum U13, torque T3 and forces (F±, F2) for an ellipse of mass mj, = 10 — 3 and moment of inertial I3 = 10 -3 using 
the old algorithm (black lines) and new algorithm (using grid Sre )• p3 an d ^2 are scaled by a factor of 100 for graphical 
purposes). The force shown on the body does not include the contribution from the external body force. 



Note that the ramp function R has three continuous derivatives since the first three derivatives of R(t) are zero at 
t = and t = 1. 

We consider an an ellipse of mass mb = 10 -3 and moment of inertia 73 = 10 -3 . The fluid is taken as an ideal gas 
with 7 = 1.4. The ellipse and fluid are initially at rest with the initial fluid state given by (p, vi,V2,p) = (I/7, 0, 0, 1). 
The smooth body force is given by ( |101[ ). The boundary conditions on the Cartesian grid, which have little influence 
for this problem, are inflow on the left with all variables given, outflow on the right side (all variables extrapolated) 
and slip walls on the top and bottom. For comparison, we solve this problem using both the old FSI algorithm and 
the new added-mass FSI algorithm. The new algorithm is run at a CFL number of 0.9. The old algorithm experiences 
exponential blowup at this CFL number and is instead run at a CFL number of 1/100. 

In the right-hand side of Fig. [l0]we show the state of the rigid body over time for the old and new algorithms. The 
body initially accelerates upward and to the right as indicated by the components of the body velocity and rotates 
in a counter-clockwise direction as indicated by the angular velocity. The forces on the body shown in Fig.[l0]do not 
include the contributions from the external body force and thus represent the force exerted by the fluid on the body. 
The force Fi indicates that the fluid pushes back on the body to nearly balance the external force ,f x (t). The results 
from the old and new algorithm are nearly indistinguishable in this plot indicating that the new algorithm provides 
an accurate approximation even with a time step that is nearly 100 times larger than the old algorithm. 

Fig. |ll| shows contours of the pressure field at times t = 0.5 and t — 1.0 for both the old and new algorithms. The 
accelerating body generates a forward moving wave that steepens over time and which has formed a shock by t — 1.0. 
The solutions from the old and new algorithm are in excellent agreement with almost no detectable differences. For 
a more quantitative evaluation of the accuracy we determine a a-posteriori error estimates by solving the problem 
on a sequence of grids of increasing resolution and using the error estimation approach described in [21, 22 . Fig. 12 



shows the estimated max-norm errors and convergence rates at time t = 0.4 when the solution is still smooth. These 
results show that the solution is converging at close to second-order accuracy. We note that for these results the 
slope-limiter was turned off in the Godunov method since this slope limiter can reduce the order of accuracy. Fig. |13| 
shows the estimated Li-norm errors and convergence rates at time t = 1.0 when the solution is no longer smooth. In 
this case the results show that the solution is converging at rates close to 1, which are the expected rates for problems 
with shocks. We note that the discrete Li-norm of a grid function is computed in the usual way by summing the 
absolute values of the values at each grid point and dividing by the total number of grid points [21] . 

8.3. Shock driven zero mass ellipse 

The shock driven ellipse problem consists of a Mach 2 shock that impacts an ellipse of zero mass and zero moment 
of inertia. This example demonstrates the robustness of the new added-mass algorithm on a difficult problem for 
which the old rigid-body FSI algorithm would fail for any time-step, no matter how small. We note that since 
the mass and moments of inertial of the body are zero in the Newton-Euler equations ( |97| l , the linear and angular 
velocities of the body respond immediately to ensure the net force on the body is zero; there is no damping in the 
response from the body's inertia. 



The overlapping grid for this problem, is the same as that used in Section 8.2 We use adaptive mesh 



refinement in some of the computations of this section. Let Qte* denote the AMR grid that has a base grid Q, 
with grid spacing As^ 3 ' « l/(10j) together with one level of refinement grids of refinement factor 4. The effective 
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Addcd-mass algorithm: t — 0.5 
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Figure 11: Accelerated ellipse: pressure at t = 0.5 and t = 1.0 for the old algorithm running at CFL number 10 2 (bottom) 
and new added-mass algorithm running at CFL number 0.9 (top) for grid Qil . 
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Figure 12: A posteriori estimated errors (max-norm) and convergence rates for the accelerated ellipse at t = 0.4 (no slope 
limiter). The scheme converges at close to second-order accuracy in the max-norm when the solution is smooth. 



resolution of the AMR grid £ r ( c X4) is thus As (jx4) « l/(40j). We note that the AMR grids are added to both the 
background grid and to the component grid around the ellipse, refer to [7] for further details of the moving-grid AMR 
approach. 

The initial conditions in the fluid consist of a shock located at x = —1 with initial state (p,u,v,p) = 
(2.6667,1.25,0,3.214256) ahead of the shock and (p,u,v,p) = (1,0,0,1.4) behind the shock. The boundary con- 
ditions are supersonic inflow (all variables specified) on the left face of the background grid and supersonic outflow 
(all vari able s extrapolated) on the other faces of the background grid. 

Fig. 14 compares the time history of the rigid body dynamics from a coarse grid, Qxt , and finer grid, C/S 2 ', 



computation. The velocity and angular velocity are seen to rapidly increase when the shock first hits the ellipse 
just after t = 0.2. The ellipse is initially accelerated up and to the right and experiences a rapid counter-clockwise 
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Figure 13: A posteriori estimated errors (Li-norm) and convergence rates for the accelerated ellipse at t = 1.0. The scheme 
converges at close to first-order accuracy in the Li-norm when the solution is not smooth. 



Shock driven ellipse, M=0 




Figure 14: Shock-drive ellipse: time histories of the center of mass, (£1,2:2), the velocity of the center of mass, (vi,V2) and the 

angular velocity 103. The colored lines are results from the coarse grid Src while the black lines are results using the finer grid 
,,(32) 



rotation. After an initial rise, the angular velocity decreases and approximately levels off at some positive valu^] 
The results from the two computations are in excellent agreement. 



Numerical schlieren and contours of the pressure field at different times are shown in Fig. 15 (see [7] for a definition 
of the numerical schlieren function). The computations were performed with AMR using the grid C/rc 6 * 4 ' (base grid 
C/re 6 ' plus one refinement level of refinement ratio 4). The solution at t — 0.4 shows the ellipse has undergone a 
rapid acceleration upward and to the right combined with a rapid counter clockwise rotation. The impact of the 
incident shock on the ellipse causes a shock to form in the region ahead of the body. By t = TO, a complex pattern of 
interacting shocks has formed in the regions above and below the ellipse. In Fig. [16] we compare the schlieren images 
of the solution at t — 1.0 from grids of different resolutions. These result show good agreement in the basic structure 
of the solution, with additional fine scale features appearing as the grid is refined. This is the expected behavior for 
inviscid computations. 



3 We note that the long time behavior of the ellipse is of interest but we do not pursue that line of investigation here. 
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Figure 15: Shock driven zero mass ellipse. Schlieren images (left column) and pressure contours (right column) at times t = 0.4, 
t = 0.6 and = 1.0 using grid qH^ x4 \ The block boundaries of the refinement grids are shown superimposed on the pressure 
contours. 
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9. Conclusions 

We have presented a stable partitioned scheme for the coupling of light rigid bodies with inviscid compressible 
fluids. This new added-mass scheme, derived from an analysis of a fluid/rigid-body Riemann problem, defines the 
force on the rigid body as a sum of the usual fluid surface forces (due to the pressure) plus an impedance weighted 
difference of the local fluid velocity and the velocity of the rigid body. The form of the added-mass terms are thus 
elucidated. The scheme was analyzed in one-dimension and shown to be well defined and stable, with a large time- 
step, even when the mass of the rigid body, ms, goes to zero. In contrast the traditional FSI coupling algorithm 
has a time-step restriction that goes to zero as nib approaches zero. Both a first-order accurate upwind scheme and 
a second-order accurate Law-Wendroff scheme were analyzed. Numerical computations in one-dimension confirmed 
the results of the theory and showed that the scheme was well behaved and accurate even when nib = 0. 

The added-mass scheme was then extended to multiple space dimensions. The result was an added-mass form of 
the Newton-Euler equations for rigid-body motion that included four added-mass tensors. The added-mass tensors 
couple the translational and angular velocities of the body and are defined in terms surface integrals involving the 
fluid impedance. Numerical results in two-dimensions were presented for both smooth and discontinuous problems. 
Second-order convergence was demonstrated using a smoothly receding piston problem with known exact solution, 
and a smoothly accelerated ellipse. The robustness of the scheme was demonstrated for the difficult case of a shock 
impacting an ellipse with zero mass and zero moment of inertia. The solution to this problem was computed on a 
sequence of grids of increasing resolution (utilizing adaptive mesh refinement), with the results on the different grids 
comparing favourably. 

There are a number of avenues open for follow-on work including the extension of the current scheme to three 
dimensions and viscous flows. In addition, we are currently investigating approaches for coupling incompressible flow 
with light bodies (both rigid and deformable). 



Appendix A. An analytic solution for the one-dimensional FSI model problem 

Consider the one-dimensional FSI problem illustrated in Fig. [2] consisting of a rigid body embedded between two 
(linearized) fluid domains. The governing equations are defined in Section [4] To simplify the presentation, we take 
the width of the rigid body to be zero, Wb = 0. The solution for Wb > follows easily from the solution for Wb = 0. Let 
the displacements in the left and right domains be defined by Ul(x, t) = J* vl(x, t) dr and Ur(x, t) = J Q * v R (x, t) dr 
respectively, and let the rigid body position be given by Ub(i). The second-order wave equations 

d tt U L (x,t)-cld xx U L (x,t) = 0, {otx<0, (A.l) 
d u U R (x,t)-c R d xx U R (x,t) = 0, forx>0, (A.2) 

describe the evolution of Ul and Ur. The evolution of the rigid body position is given by the rigid body equations 
of motion with the applied stress from the fluid determining the force on the body 

m b dttU b (t) = p R c 2 R d x U R (0,t) - p L cld x Ui(0,i). (A.3) 

Assume given initial conditions 

U L (x,0) = U (x), d t U L (x,0) = Vo(x), fora;<0 (A.4) 

U R (x,0) = U (x), d t U R (x,0) = Vo(x), foTx>0 (A.5) 

U b {0) = U (0), 9*14(0) = Vo(0). (A.6) 

The exact solution for x < can be written in terms of the d'Alembert solution as 

U L (x,t) = f L (x- c L t) + g L (x + c L t), (A.7) 

where 

h(0 =i (u (0 - - ^ V (s) ds) , (A.8) 



2 V *. J 

»<o-{ + £ J? *<•>*)■ f »«°' (A , 9) 

1 fi(4)-A(-{). far « <0. 

Likewise for x > 0, the solution can be written 

U R (x, t) = f R (x - cat) +g R (x + c R t) (A.10) 
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where 




for £ > 
for £ < 0, 



U (0 + 



cr 



V (s)ds 



(A.ll) 
(A.12) 



For x < —CLt or x > crI the solution is given by the usual d'Alembert solution for the Cauchy problem, 

fX-\-ct 



U(x,t) = ^(Uo(x-ct) + U (x + ct)) + — 



Vo(s) ds, 



where c = cl or c = cr for the left and right domains, respectively. For —erf < x < crX, the left and right solutions 
are coupled to the rigid body. For this case, the unknown interface position Ub is found as the solution to the linear 
ODE 

m b duU b (t) + {z R + z L ) d t U b (t) = g{t) (A.13) 

where g(t) — pRC 2 R d x Uo(c R t) — pLC 2 L d x lIo(—c L t) + z R Vo(c R t) + z L Vo(—c L t). Solutions to the corresponding homoge- 
neous ODE mbdttrj(t) + (zr + zl) dti](t) = are easily found as 

r ll (t)=e' t(ZR+ZL)lm \ and ^(i) = 1. 

The method of variation of parameters can be used to derive an exact solution to ( |A.13[ ) by looking for a solution of 
the form 

U b (t) = fci(*)?7i(t) + k 2 {t) V2 (t). (A.14) 
The unknown functions k± (t) and k 2 (t) are found to be 



k 2 (t) 



m(t)g{t) 

W[ Vl , V2 ](t) 

m{t)g{t) 



dt + const, 



w[rn,m](t) 



dt + const, 



(A.15) 
(A.16) 



where rj 2 ]{t) is the Wronskian of the homogeneous solutions. The integration constants are determined by the 

initial conditions. For a more detailed discussion on solution methods for ( |A.13[ | refer to for example. 

A specific solution of the form (A.14I is determined by specifying initial conditions Uo(x) and Vo(x). We illustrate 
with an example where an initial Gaussian pulse (of velocity and stress) moves from left to right and interacts with 
the rigid body and fluid domains as time progresses. Let the initial conditions be given as 



U (x) 
V (x) 



1 V^rerf (/3(x - x )) 
4~ 



CL 



exp (— ft {x 



- x 



(A.17) 
(A.18) 



Here /3 > and xq < are parameters used to define the center and width of the initial pulse. Also notice that we 
envision the pulse to originate entirely in the left domain which is the reason for the appearance of cl in the initial 
condition definition. The velocity of the rigid body can be found as 



U b (t) = 



zr{cr - c L )^/ n 
kcRfirnb 

erf 



■ exp 



(z L + Zr)(z L + Zr- 4/3 2 m b CR(CRt - x )) 



4c 2 R m b 2 /3 2 



zl + z R — 2cRf3 2 nib(cRt — x ) 



erf 



2/37716 



exp 



2cRirib/3 

(z L + z R )(z L + zr — ic L /3 2 m b (cLt + x'o)) 
Ac 2 L m b 2 P 2 
zl + zr — 2cLP 2 m b (c L t + x 



CL 



exp 



erf 



Q 2 2 

-P x 



(zl 



2c L m b l3 
ZR)t s 



erf 



zl + zr + 2cR/3 2 mbXo 
2c R m b P 



ZL + ZR — 2clP 7776X0 



2c L m b 



nib 



(A.19) 



Analytic expressions for the position and acceleration of the body are determined from (A.19 1, by integration and 
differentiation, respectively. Note that (A.19 1 is not easily evaluated numerically with standard math libraries as 
777,;, — > 0. For the small mass case, ( |A.19[ ) can be evaluated using asymptotic expansions of the error functions as their 
arguments approach plus or minus infinity. The desired level of accuracy can be obtained by appropriately truncating 
the resulting series expansion. In practice, we find that for 7776 ~ 0.1 such a procedure should be used. 
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Appendix B. Examples of added mass matrices for constant fluid impedance 



In this section we illustrate the form of the added mass matrices defined by (93 1- (94 1, for some common body 
shapes when the fluid impedance 2/ is taken to be constant. We denote the entries of A vv by af/, the entries of A vt * 
by a%f and the entries of A uu by af". 

Appendix B.l. Added-mass matrices for an ellipse 

Consider a two dimensional ellipse with semi-axes of length a and b and center of mass xo = 0. A point on the 
ellipse is x(#) = [acos(8), bsin(8), 0] T . The tangent to this point is 



Xfl 

ll x s|l 



-asin(0), 6cos(0), 0] T / J a 2 sin 2 (6) + b 2 cos 2 (8). 



Thus 



n= [6cos(0), asin(6>), 0] T / \J a 2 sin 2 (8) + b 2 cos 2 (8), 
y— [acos(8), bsinS, 0] T , 



and 



Yn=[0, 0, (a 2 - b 2 ) cos(6>) sin(0)] T /yV sin 2 (6>) + b 2 cos 2 (6>). 
Thus (leaving out some zero rows and columns which do not apply in two-dimensions) 

T 1 



Yn(Yn) J = 



a 2 sin 2 (0) + 6 2 cos 2 (8) 



a 2 sm 2 (8) + b 2 cos 2 (8) 



6 2 cos 2 (0) a6cos(l9)sin(6>) 



abcos(8) sin(0) 



■(0) 






(a 2 - & 2 ) 2 cos 2 (60sin 2 (6>) 



The increment in arclength is ds = Vdx ■ dx = y/ a 2 sin 2 (8) + b 2 cos 2 (8) dd. Thus 



Jo \fa- 



2 / 



y/a 2 sin 2 (8) + b 2 cos 2 (8) 



6 2 cos 2 (0) afocos(0)sin(0) 
a&cos(0)sin(0) a 2 sin 2 (0) 



dO. 



(B.l) 
(B.2) 

(B.3) 



A" w = 



and 



A vu = {A uv y = 



























033 




"0 






















^a 2 sin 2 (6») + b 2 cos 2 (8) 






(a 2 - & 2 ) 2 cos 2 (0)sin 2 (0) 



dO. 







z f 

y/a 2 sin 2 {8) + b 2 cos 2 (8) 



b(a 2 - b 2 )cos 2 (0)sin(0) 
a(a 2 - 6 2 )cos(0) sin 2 (8) 




dO. 



(B.4) 



(B.5) 



Values for a™, a^, and a^, (which can be written in terms of elliptic integrals) for some ratios of b to a are given 
in Figure B.17 The values for a^li a i3 an d 0-23 are zero f° r uniform Zf (but can be non-zero when zj varies). Note 
that for the case of a circle, a — b, aii = 022 = (zf /a)ira 2 where ira 2 is the area of the circle. Compare this result to 
that for the sphere in Section [Appendix B.2| 





b = a 


6 = a/2 


6 = a/10 


6 = a/100 


vv 

"11 


TTZfa 


1.26z/ a 


.108.3/0 


.0020z/a 


a 22 


TTZfa 


3.58z/a 


3.96z/a 


3.99z/a 


"33 





.581z /a 3 


l^z/a 3 


1.330 / a a 



Figure B.17: Components of the added-mass matrices for an ellipse for various values of b/a with constant Zf. Values for 
b/a^l are approximate. 
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Appendix B.2. Added-mass matrices for an ellipsoid 

We consider an ellipsoid with semi-axes of length a, b and c and center of mass at xo = 0. A point on the surface 
of the ellipsoid is given by 



x.(9,4>) = [asm(<p) cos(#), 6sin(</>) sin(0), ccos(0)] 



€[0,tt], 6>€[0,2tt]. 



From this formula it is straightforward to determine n and Yn in the formulae for the added mass matrices. For a 
sphere of radius a, i.e. a — b — c, we get (4tt/3 « 4.18879) 



.4" 



4tt/3 
4tt/3 
4tt/3 



z f a 









z f a 









(B.6) 



Recall that the volume of the sphere is V — 47ra 3 /3 so that a™ = (zf/a)V, i = 1,2,3. The rotational added-mass 
entries af" , i = 1,2,3 are zero since a rotating sphere exerts no force on the adjacent (inviscid) fluid. 

For b = a, c = 2a, we can compute the added-mass matrix entries approximately by quadrature giving the values 



z t a 



9.254 
9.254 
2.971 



A" 



z f a 









A u 



z f a 



4.712 
4.712 




(B.7) 



This ellipsoid is longest along the z-axis and has circular cross-sections for z constant. The values of dn and 022 
are larger than 033 which indicates that the added mass is larger for translational motions in the x- or y-directions 
compared to the z-direction. In other words is takes more force to move the ellipsoid in the x— or y— directions 
compared to the z-direction. This is consistent with the shape of the ellipsoid which is longest along the z-axis and 
thus has a greater effective cross-sectional area when viewed from the x— or y— directions. 
For b = 2a, c = 3a, the added-mass matrix entries are given approximately by 



* vv 2 

A = Zfa 



32.307 
11.023 
5.552 



A vu = Zfa 3 - 









A^ = z f a 4 



6.840 







53.511 







15.963 



(B.8) 



In this case, the translational added mass a™ is largest, consistent with the effective cross-sectional area being largest 
when the ellipsoid is viewed in the :r-direction. In other words it takes more force to move the ellipsoid in the 
a;-direction, compared to the other directions. 



Appendix B.3. Added-mass matrices for a rectangle 

The added-mass matrices for bodies with piecewise constant surface normals are easily computed. Consider the 
rectangular body of length l x , height l y , and center of mass x = given by 1Z = {(x, y) \ — l x /2<x< l x /2, — l y /2 < 
V < ly/2}- The added-mass matrices for this case are 



A vv = z f 



2ly 







2l x 




A™= zf 









A™ = z f 






i(Jj»+J 



(B.9) 



Appendix B.4- Added-mass matrices for a rectangular prism 

Finally, consider the rectangular prism with dimensions l x , l y 
{(x, y,z) \ — l x /2 < x < l x /2, — ly/2 < y < ly/2, — l z /2 < z < l z /2, }. The added-mass matrices are 



and center of mass x = given by V = 



z f 



2Zy I z 










X>>Z 





-I .! 



A v 



z f 



A u 





(l 3 x + It) 




l t(l 






3 +^) 



(B.10) 
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